Free energy determination of phase coexistence in model Ceo: 
A comprehensive Monte Carlo study 
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The free energy of the soUd and fluid phases of the Girifalco Ceo model are determined through 
extensive Monte Carlo simulations. In this model the molecules interact through a spherical pair 
potential, characterized by a narrow and attractive well, adjacent to a harshly repulsive core. We 
have used the Widom test particle method and a mapping from an Einstein crystal, in order to 
estimate the absolute free energy in the fluid and solid phases, respectively; we have then deter- 
mined the free energy along several isotherms, and the whole phase diagram, by means of standard 
thermodynamic integrations. The dependence of the simulation's results on the size of the sample 
is also monitored in a number of cases. 

We highlight how the interplay between the liquid-vapor and the liquid-solid coexistence con- 
ditions determines the existence of a narrow liquid pocket in the phase diagram, whose stability 
is assessed and conflrmed in agreement with previous studies. In particular, the critical tempera- 
ture follows closely an extended corresponding-states rule recently outlined by Noro and Frenkel [J. 
Chem. Phys. 113, 2941 (2000)]. 

We discuss the emerging "energetic" properties of the system, which drive the phase behavior in 
systems interacting through short-range forces [A. A. Louis, Phil. Trans. R. Soc. A 359, 939 (2001)], 
in order to explain the discrepancy between the predictions of several structural indicators and the 
results of full free energy calculations, to locate the fluid phase boundaries. 

More generally, we aim to provide extended reference data for calculations of the free energy of 
the Ceo fuUerite in the low temperature regime, as for the determination of the phase diagram of 
higher order Cn>60 fuUerenes and other fullerene-related materials, whose description is based on 
the same model adopted in this work. 
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PACS numbers: 61.48. -fc, 64.70.Fx 



I. INTRODUCTION 



We report an extensive investigation of the free energy 
characteristics of the Girifalco model of Cgo fullerene [1] . 
As is well known, this representation hinges on the fact 
that Geo molecules have almost spherical shape and freely 
rotate at sufficiently high temperatures [2]. Under these 
conditions, the hollow molecular cages can be assimilated 
to spheres whose surface consists of a uniform distribu- 
tion of carbon sites. The overall interaction between two 
fullerene particles is then obtained by an integral of the 
interaction between pairs of sites on different cages, even- 
tually yielding an analytic central two-body potential. 
The latter is characterized by a harshly repulsive core at 
short distance, followed by a deep attractive well which 
rapidly decays with the interparticle distance [1]. 

The Girifalco model constitutes a prototype system 
in several respects; recent studies suggest that a similar 
"smeared out" spherical description can be attempted for 
Cn>60 systems, although fullerene molecules with n > 60 
can have a sensibly non-spherical shape; the cases n = 70, 
76 and 84 have in particular been examined (see [3] and 



references therein). On the other hand, more refined cal- 
culations of the fullerene-fullerene interaction yield re- 
sults very close to those predicted through the Girifalco 
model [4,5], and a similar representation has been used 
for the description of other hollow nanoparticles as car- 
bon onions, or metal dichalcogenides (also termed inor- 
ganic fuUerenes) as GaAs and CdSe [6]. Moreover, a 
modification of the Girifalco model, suitable for the de- 
scription of solid Ceo at low temperatures, has been re- 
cently proposed [7] ; this development seems of particular 
interest since fuUerites, doped with organic molecules and 
upon the injection of electron (or holes), exhibit a super- 
conducting behavior up to T = 112 K [8]; an accurate 
description of such a simple model might prove useful 
for further studies on the lattice behavior upon impurity 
doping. 

In our opinion, such a possible reference role of the Gir- 
ifalco model for further studies on fuUerenes and other 
systems, calls for a complete and confident determina- 
tion of its phase diagram. With this purpose in mind, we 
have investigated the free energy characteristics of the 
Geo model for both the solid and the fluid phase through 
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extensive Monte Carlo simulations, spanning the whole 
high-temperature region of the phase diagram. 

A second general aspect of the Girifalco potential is 
related to its short-range nature. It is useful to recall 
in this respect that a ubiquitous definition of the range 
of interaction has been recently proposed in Ref. [9] . As 
several studies have pointed out, the most apparent con- 
sequence of a reduced interaction length is the metasta- 
bility of the liquid- vapor equilibrium, which is preempted 
by the fluid-solid coexistence [10]. It is actually known 
that a stable liquid-vapor coexistence still survives for 
the model envisaged here, albeit restricted to a few tens 
degree temperature range [11,12]. A tiny liquid pocket 
has been predicted also in the phase diagram of similar 
models of higher order fuUerenes [3], although it is ar- 
gued that even a modest reduction of the range of the 
forces might cause the disappearance of the liquid phase. 
The emerging borderline nature of the Girifalco model 
implies that the overall appearance of its phase diagram 
sensitively depends on the details of the interaction po- 
tential; in fact, the initial controversy around the exis- 
tence of a stable liquid phase for this system [11,13] has 
been solved by taking into account on one side the full 
role of the attractive part of the interaction [14], and on 
the other, a conveniently large simulation sample [15]. 
More generally, it has been recognized in Ref. [16] that 
the physical behavior of a wide class of systems charac- 
terized by short-range interactions, and in particular the 
onset of freezing, are substantially affected by the "per- 
turbative" part (with respect to the repulsive core) of 
the potential, rather than being dominated by cxcluded- 
volume and packing (i.e. by entropic) effects, as is the 
case in the currently accepted van der Waals picture of 
simple liquids [17]. We here investigate, as a further key 
purpose of this work, the same issue for the Cgo model 
and try to reconcile the predictions coming from differ- 
ent structural indicators usually related to th(^ frc(ving 
threshold of several simple liquids [18,19] and adopted in 
early calculations for this model [11,20,21] — with the 
results of full free energy calculations. 

The paper is organized as follows: we present in Sec- 
tion II the model and the simulation strategics: the 
results are reported and discussed in Sect. Ill, while 
Sect. IV is devoted to the conclusions and future per- 
spectives of our investigation. A preliminary account of 
this work has been presented elsewhere [22]. 



II. MODEL AND SIMULATION STRATEGIES 

The Girifalco potential between two Ceo molecules, is 
written as [1]: 



v{r) = —a\ 
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where s = r/d, ai = N'^A/I2d^, and = N"^ B / ^M'^'^ : 
N = 60 and d = 0.71 nm are the number of carbon atoms 
and the diameter, respectively, of the spherical particles; 
A = 32x 10-60 ergcm^ and B = 55.77 x 10-1°^ ergcm^^ 
are constants entering the Lennard- Jones 12-6 potential 
through which two carbon sites on different molecules are 
assumed to interact. The finite distance at which the po- 
tential (1) crosses zero and the minimum of the potential 
well depth arc a ~ 0.959 nm and e ~ 0.444 x 10"-^^ erg 
at Tmin = 1-005 nm, respectively. 

In order to investigate the coexistence properties of 
model (1), it is required the knowledge of the free ener- 
gies of both the solid and the fluid phase. As a general 
strategy, the free energy of the system is first evaluated 
all along a supercritical isotherm at temperature T by 
integrating the pressure P as a function of the density p 
according to the formula [23] : 
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here FjlSS is the Helmholtz free energy per particle, /3 = 

l/ksT is the inverse temperature, ks is the Boltzmann 
constant and (p, T) is a thermodynamic state where the 
free energy is known (see below). The free energy at 
different temperatures is then calculated along isochoric 
paths as: 



l3F{p,T) _ PF{p,T) 



N 
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(3) 



where U/N is the internal energy per particle of the sys- 
tem. The equilibrium conditions finally derive from the 
equality of the chemical potential and the pressure of the 
different phases. 

All quantities entering Eqs. (2) and (3) have been 
determined through standard Monte Carlo simulations 
at constant volume or pressure. Simulations have been 
mostly carried out on a sample composed of N = 864 Ceo 
particles enclosed in a cubic box with periodic boundary 
conditions. The Ceo interaction has been considered up 
to half the box length. 

As far as the determination of the reference free en- 
ergy F{p,T) in Eq. (2) is concerned, we have used for 
the solid phase the Einstein crystal method described by 
Prenkel and Ladd [24,25]. Namely, the Cgo interaction 
in the solid fullerito is smoothly transformed, through 
a coupling parameter A, into a corresponding harmonic 
potential so that the configurational energy U ({r}) takes 
the following form [25]: 

U{M) = UcoiiM) + (1 - A) [Uc,o{{r}) - f/ceo({ro})] 



N 



(4) 
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where ro,i is the lattice position of atom i and t/ceo ({'"o}) 
is the static contribution to the potential energy; a is the 
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spring constant of the Einstein crystal. The free energy 
difference can then be written as; 

pF{-p,T) _ /3Fceo 
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where the configurational free energy of the Einstein crys- 
tal is: 
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In the equation above, the last two terms account for the 
fixed center of mass constraint at which simulations are 
performed [25]. 

As for the absolute free energy of the fluid phase, the 
chemical potential jj, has been estimated at several in- 
termediate densities through the Widom test particle 
method [25,26]. 

In parallel with the free energy, the entropy of the fluid 
phase is also systematically analyzed, given that in sev- 
eral earlier papers the onset of freezing of the Girifalco 
model has been associated with the vanishing of the resid- 
ual multiparticle entropy As, namely 



As = Sc 



S2=0, 



(7) 



according to the one-phase freezing criterion originally 
proposed in Ref. [19]. In Eq. (7), Sex is the the excess 

entropy per particle of the system (in fcs units) and S2 is 
defined in terms of the radial distribution function g{r) 
of the system as [27]: 



S2 



-| J{g{r)\n[g{r)]-g{r) + l}dr. 



(8) 



In previous works As has been evaluated using liquid 
state integral equation theories. Here we report Monte 
Carlo data for the structural and thermodynamics quan- 
tities entering Eqs. (7) and (8), in order to perform a 
rigorous test of the theoretical predictions based on the 
criterion (7). 



III. RESULTS AND DISCUSSION 

The solid and liquid branches of the equation of state, 
to be integrated in Eq. (2), have been calculated at the 
supercritical isotherm T = 2100 K through both NVT 
and NPT Monte Carlo simulations, in order to cross- 
check the predictions of the two algorithms. Five to six 
runs of 25 000 steps have been performed at each thermo- 
dynamic point investigated. The compressibility factor 
(3P/p is shown in Fig 1. 



Results for the absolute free energy of the solid phase 
at T = 2100 K and p = 1.375 nm^'', as obtained through 
the Einstein crystal method, are reported in Fig. 2. Con- 
stant volume simulations have been carried out for sev- 
eral values of the switching parameter A. The spring 
constant is set to a value a/e = 490, which makes 
the interactions in the pure Einstein crystal as close 
as possible to those of the original system, so to opt- 
mize the accuracy of the numerical integration scheme 
of Eq. (5) [25]. We have analyzed the free energy de- 
pendence on the system size; results with 256 and 2916 
particles are shown in Fig. 2. It appears that the ef- 
fect of N on the estimate of the free energy is small but 
systematic; in the thermodynamic limit {N cx)) the 
smooth extrapolation in the bottom panel of Fig. 2 yields 
(3F/N{p = 1.375 nm-3,T = 2100 K) ~ -1.401. 

As for the fluid phase, the excess chemical poten- 
tial at T = 2100 K and p = 0.60 nm'^ has been 
calculated through the Widom test particle method as 
/3/iox = —1.523. Several tests at higher densities have 
also been conducted, in order to assess the results ob- 
tained via Eq. (2). The free energy at T = 2100 K is 
shown in the top panel of Fig. 3, along with the common 
tangent construction, which determines the coexistence 
conditions. The fi vs P behavior is displayed in the bot- 
tom panel of Fig. 3; it emerges that the thermodynamic 
integration of Eq. (2) fully agrees with the direct estimate 
of the chemical potential based on the Widom technique. 

Starting from the knowledge of the free energy along 
the isotherm T = 2100 K, the free energy at differ- 
ent temperatures has been obtained through Eq. (3). 
We have examined the isochores p = 1.25, 1.27, and 
1.30 nm~^ in the solid phase and the density range 
p = [0.70 — 1.00] nm~^ with steps Ap = 0.05 nm~^ in 
the liquid phase, descending down to T = 1800 K, with 
temperature intervals AT = 25 K. Simulations have been 
carried out with N = 864 particles at constant density; 
four to eight cumulation runs of 10 000 steps at each state 
point are sufficient to yield accurate internal energies es- 
timates (see Fig. 4), so to allow a smooth interpolation 
for the integration in Eq. (3). The pressure and the free 
energy along several isotherms are shown in Figs. 5 and 6, 
respectively; fully consistent results are obtained if we de- 
termine the free energy along an isochoric path first, and 
then integrating Eq. (2) at constant temperature. The 
accuracy of this global check is also evidenced in Fig. 6. 

In order to analyze the thermodynamic properties 
of the system in the low density regime p = [0.05 — 
0.20] nm~^, we have calculated the chemical potential 
of the vapor phase along the isotherm T = 1900 K 
through the Widom technique, and then estimated 
through Eq. (3) the free energy in the temperature range 
T = [1800 - 1900] K. 

As is visible in Fig. 7, where the chemical potential 
of the different phases is shown, liquid-solid equilibrium 
is stable at T = 1900, while a solid-vapor coexistence 
takes place at T = 1850 K. The intermediate tempera- 
ture T = 1875 K is characterized by almost a compara- 
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ble value of the chemical potential of the various phases. 
We thus estimate that the triple point temperature is 
Ttr ~ 1875 K, at the pressure Ptr ^ 2.4 MPa; we then 
obtain from the equation of state ptr — 0.74 nm^'', in 
close agreement with the free energy results of Ref. [12], 
T = 1880 K and p = 0.74 nrn'^. The features in the /j. 
vs P behavior which determine the narrow temperature 
width of the liquid pocket are clearly illustrated in Fig. 7, 
where it emerges that the solid and the liquid free energy 
branches have dissimilar slopes and considerably differ- 
ent spacing under equal variations of temperature. These 
two circumstances cause a fairly rapid shift of the liquid- 
solid intersection points and hence of the corresponding 
equilibrium parameters. Conversely, the vapor branch is 
hardly sensitive to temperature variations, and already at 
T > 1900 K tends to loose any further intersection with 
both the solid and the liquid branch. The only surviving 
liquid-solid equilibrium therefore fully characterizes the 
system behavior for T > 1950 K. 

The phase diagram of the system is displayed in Fig. 8, 
along with our previous Gibbs Ensemble Monte Carlo 
(GEMC) determination of the binodal line [15]. Dis- 
tinct liquid-vapor and liquid-solid equilibria take place 
at temperatures slightly higher than 1875 K; as far as 
the liquid-vapor coexistence is concerned, we observe a 
remarkable agreement between free energy calculations 
and GEMC results. At lower temperatures, the binodal 
points are metastable with respect to the freezing line; 
in this well, the GEMC approach fully reproduce 

the free energy data. We thus retain the GEMC results 
T„. ~ 1940 K, Per 0.42 nm~3, and Per ^ 2.7 MPa, 
obtained in Ref. [15] with 1500 Ceo particles, as a reli- 
able estimate of the critical point parameters. We note in 
this context that an extended correponding-state behav- 
ior, which states a linear relationship between the critical 
temperature and the range R of the interaction potential, 
has been recently outlined in Ref. [9]. R is defined by a 
map of the potential into an effective square well inter- 
action with the same second virial coefficient. We obtain 
for the Ceo model R = 0.16, a value immediately over 
the minimum threshold of a stable liquid- vapor equilib- 
rium, R = 0.13 — 0.15 [9]; the extended rule's estimate 
for the critical temperature is then T = 1922 K, in fair 
agreement with the above Tci = 1940 K result. 

A comparison with the phase diagram determined by 
Hagen and coworkers [13] is reported in Fig 8. It has 
been conjectured, on the basis of a theoretical density- 
functional investigation [14], that the emerging discrep- 
ancy between Ref. [13] results and subsequent calcula- 
tions might be due to an early cutoff of the Cgo in- 
teractions, which substantially affects the location of 
the liquid-vapor binodal line. We argued on the other 
hand [15], that a fairly large simulation sample should 
be employed in this case, in order to take into account 
the peculiar density fluctuations in the GEMC simula- 
tions, which act to destabilize the liquid- vapor separa- 
tion. We note that, while the overall fluid phase bound- 
aries can sensitively depend on the interaction details 



(and hence on the cutoff), truncation effects play a minor 
role in the determination of the melting line; the latter 
appears indeed coincident with our estimate, and almost 
independent on temperature variations, ranging from 
p = 1.26 nm-3 at T = 2200 K to p = 1.28 nm'^ at T = 
1800 K. Our results positively agree with the phase dia- 
gram obtained in Ref. [12] (also shown in Fig. 8), where 
a slightly higher critical temperature T = 1954 — 1980 K 
is reported. 

The phase diagram in the P — T representation is dis- 
played in Fig. 9, where it appears that also the pres- 
sure range of the liquid phase is rather restricted, span- 
ning only a few tens bar over the triple point pressure 
Ptr = 2.4 MPa. 

We now turn to the indications on the freezing condi- 
tions obtained through the one-phase criterion expressed 
by Eq. (7). As is visible in Fig. 8, the As = locus tends 
to overestimate the coexisting fluid density, thus affecting 
the location of the triple point in the phase diagram. We 
remark that a similar trend would emerge if the Hanscn- 
Verlet prescription [18,23] for the height of the first peak 
of the structure factor were used [3]. In fact, the close 
correspondence between the two interpretations of the 
freezing transition in terms of structural indicators, has 
been recently pointed out in Ref. [28]. Similarly, the be- 
havior of the internal energy, as well as the height of 
the first peak of the radial distribution function, show a 
non-monotonic behavior [21], suggesting that the system 
becomes unstable against the phase separation around 
the density p ~ 1.0 nm~^ (see also Ref. [29]). 

It thus appears that such indicators identify a very 
restricted range (if not a unique locus) of density vs tem- 
perature states, over which a structural reorganization 
of the fluid phase should be tendentially established, in 
order to satisfy purely excluded- volume, or equivalently 
entropic, demands. The almost vertical disposition of 
the As = line (see Fig. 8), which means that the den- 
sity keeps constant along the locus irrespective of the 
temperature, clearly reflects the substantial absence of 
any energy scale associated with such an indication. It 
is known at present that the vanishing of the residual 
multiparticle entropy accurately predicts the thermody- 
namic freezing threshold for a wide class of simple fluids, 
including hard-core models and the Lennard-Jones po- 
tential, in both three and two dimensions [28]. In these 
models, the interparticle potential is dominated by steric 
effects, while the attractive forces can be treated as a 
perturbation to the inherent hard-sphere system, which 
essentially drives the liquid behavior, an approach com- 
monly called the van der Waals picture of fluids [17]. 
This is furthermore illustrated in Fig. 8, where it emerges 
that the reference hard-sphere system for the Ceo model, 
obtained by splitting the potential into a repulsive and 
a perturbative part in the WCA fashion (after Weeks et 
al [30]), also freezes around the locus of vanishing residual 
multiparticle entropy, p(T) ~ 1.0 nm~^. 

We argue that the discrepancy between one-phase indi- 
cators and full free energy calculations about the freezing 
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transition of the system at issue, leads to a different sce- 
nario, where the strong attractive and rapidly decaying 
well in the interaction potential critically affects the con- 
ditions for the onset of the solid-liquid transition. The 
latter is driven in this case by "energetic" rather than 
by entropic effects, in agreement with several indications 
recently collected for systems with short-range interac- 
tions by A. A. Louis [16]. It is shown in Ref. [16] that 
for this class of fluids most of the features exhibited by 
hard-sphere dominated systems do not arise, resulting in 
particular in the anticipation of the frezing threshold to 
lower densities than those predicted by solely structural 
conditions. In this respect, the lack of an accurate esti- 
mate of the freezing line of model Ceo upon use of either 
the residual multiparticle entropy or the Hansen- Verlet 
prescriptions, must be interpreted less a "failure" of the 
criteria themselves, than a manifestation of their limited 
applicability in the present context, both indicators pre- 
dicting the freezing threshold of the fluid on the basis of 
almost purely entropic requirements. 

It appears in conclusion that the phase behavior of 
the Girifalco Ceo model can be consistently understood 
in terms of general properties of systems interacting 
through short-range forces. Nevertheless, several issues 
are still open to further investigations and we refer firstly 
to the physical meaning to be associated to the fluid 
phase boundaries signalled by the above structural in- 
dicators. We note, on the other hand, the lack of clear 
indications of the freezing transition at a thermodynamic 
level, both in the internal energy and in the pressure, as 
documented in Figs. 4 and 5. More generally, the ques- 
tion of a detailed description of the microscopic behavior 
of the Ceo model — as of other fluids interacting through 
short-range forces already raised in Ref [31], is still un- 
solved. The formation of metastable clusters of strongly 
correlated particles in such fluids has been discussed in 
some recent papers [32] , on the basis of the short-range 
nature of the interaction potential; however, in a recent 
molecular dynamics study [33] of quite a similar model, 
we have not been able to identify any net precursor, at 
a microscopic level, either of the fluid-solid threshold, or 
of the incipient crystallization of the liquid phase. 



IV. CONCLUSIONS AND PERSPECTIVES 

The free energy of the solid, liquid and vapor phases 
of the Girifalco model of Ceo has been studied by means 
of extensive Monte Carlo simulations at constant den- 
sity or pressure; the full phase diagram of the system 
has then been reconstructed on such a basis. It is con- 
firmed by this comprehensive investigation that a stable 
liquid phase for this system exists, albeit confined to a 
rather restricted temperature range: we confidently esti- 
mate the triple and critical temperatures as Ttr — 1875 K 
and Tcr = 1940 K, respectively. The pressure range of the 
liquid phase also appears rather narrow, spanning only 



a few tens bar over the triple point pressure. Our re- 
sults illustrate how the interplay of various free energy 
branches determines the overall appearance of the phase 
diagram, and in particular the narrow extension of the 
liquid pocket. 

The estimate of the freezing conditions, based on sev- 
eral structural indicators, is also critically discussed. 
It turns out that such indicators identify a practically 
unique thermodynamic locus, where the fluid phase be- 
comes unfavoured due to entropic requirements. This lo- 
cus almost coincides with the true thermodynamic freez- 
ing line for systems whose phase behavior is dominated 
by steric effects. For the model at issue, the solid-liquid 
transition is instead strongly affected by the deep, short- 
range attractive well in the interaction potential. As a 
result, the freezing transition of the fluid is driven to 
lower densities, mainly by energetic effects, in agreement 
with a scenario recently proposed by other authors. 

As far as further studies are concerned, we are cur- 
rently assessing, against the wide set of data produced 
in this work, the performances of several refined integral 
equation theories of the liquid state [34] and perturba- 
tion approaches, in order to describe on a full theoretical 
ground the phase diagram and the free energy properties 
of the Girifalco model. A refinement of our preliminary 
report on the phase diagram of higher order fuUerenes 
Cn>60 [3] is also in progress. Finally, we plan to inves- 
tigate the Ceo properties, as well as the phase behavior 
of systems constituted by doped fullerites, in the low- 
temperature region of the solid phase. 
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FIG. 1. Equation of state of the Girifalco model in the fluid (circles) and solid (squares) phases at T = 2100 K, as obtained 
through NVT (solid symbols) and NPT (open symbols) Monte Carlo simulations. 
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FIG. 2. Results of the Einstein crystal procedure for the determination of the free energy in the solid phase at T = 2100 K 
and p = 1.375 nm~^. Left: integrand in Eq. (5) (circles), resolved into the harmonic (downward triangles) and Ceo (upward 
triangles) contributions, see Eq. (4). Simulation results obtained with N = 256 (diamonds) and N = 2916 (crosses) particles 
arc also shown. Lines arc smooth interpolations of the data points. Right: size dependence of the free energy of the Ceo crystal; 
the diamond indicates the extrapolation (solid line) to the thermodynamic limit. 
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FIG. 3. Free energy (left) and chemical potential (right) of the Ceo model along the isotherm T = 2100 K. Squares: solid 
phase; circles: fluid phase. Solid lines are smooth interpolations of the data points. In the left panel the common tangent 
construction (dashed line) is shown. In the right panel the diamond locates the coexistence conditions; the direct estimates of 
the chemical potential (triangles), based on the Widom test particle method, are also reported. 
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FIG. 4. Internal energy per particle in the fluid (top) and solid (bottom) phases. The behavior along the isochores p = 0.80, 
0.85, 0.90, and 0.95 nm~^ (top panel, from top to bottom) and p = 1.251, 1.27, and 1.305 nm~^ (bottom panel, from top to 
bottom) is shown. Lines are smooth interpolations of the data points. 
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FIG. 5. Equation of state in the fluid (left) and solid (right) phases. The behavior along the isotherms T = 2200, 2100, 2000, 
1950, and 1900 K (left panel, from top to bottom) and T = 2200, 2100, 2000, 1950, 1900, 1850, and 1800 K (right panel, from 
top to bottom) is shown. Lines are smooth interpolations of the data points. 



ea 




-2.6 




0.76 



0.78 



0.8 1 



1.1 



1.2 



1.3 



1/p [nm ] 



FIG. 6. Helmholtz free energy per particle in the solid (left) and fluid (right) phases. The behavior along the isotherms 
T = 2200, 2150, 2100, 2050, 2000, 1950, 1900, and 1875 K (from top to bottom) is shown. Solid hnes are guides to the eye; 
dashed lines are obtained by integrating Eq. (2), see text. 
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FIG. 7. Chemical potential vs pressure in proximity of the triple point, in the vapor (V, full lines), liquid (L, dot-dashed lines) 

and solid (S, dashed lines) phases at T = 1900, 1875, and 1850 K. Liquid-solid (squares), liquid- vapor (circles) and solid-vapor 
(triangles) coexistence points are shown; full and open symbols refer to stable and metastable equilibria, respectively. 




FIG. 8. Phase diagram of the Girifalco model according to the free energy investigation of this work (solid squares, coexistence 
points; open squares, metastable liquid- vapor separation). The line represents the Gibbs Ensemble Monte Carlo predictions for 
the liquid- vapor coexistence [15]. The critical (full diamond) and triple (open diamond) points are also shown. Open upward 
triangles and circles are the simulation results of Ref. [12] and [13], respectively. Open downward triangles: As = locus; solid 
downward triangles: freezing line of the hard-sphere fluid corresponding to the Ceo model through the WCA prescription (see 
text). 
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FIG. 9. Phase diagram in the P — T representation. Solid squares, this work; open squares, Gibbs Ensemble liquid-vapor 
coexistence points [15]. The triple (open diamond) and critical (solid diamond) points are also shown. Lines are intended as 
guides to the eye. 
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